home *** CD-ROM | disk | FTP | other *** search
/ Ham Radio 2000 #2 / Ham Radio 2000 - Volume 2.iso / HAMV2 / ANTENNA / YAGIU112 / GAIN.C < prev    next >
C/C++ Source or Header  |  1995-08-07  |  4KB  |  132 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include <malloc.h>
  4. #include <errno.h>
  5. #include "yagi.h"
  6. #define TINY 1e-10
  7.  
  8. extern int errno;
  9.  
  10. /* This function finds the gain both in the E plane (xz plane) and the
  11. H plane (xy plane) at angle (theta, phi). The method used is as described
  12. on page 1-12 of 'Yagi Antenna Design' by Dr. Lawson , ARRL */
  13.  
  14. void gain(double theta, double phi, double pin, double F, struct element_data *coordinates, struct FCOMPLEX *current, int elements, double *gain_E_plane, double *gain_H_plane, double actual_frequency, double design_frequency) 
  15. {
  16.     int i;
  17.     double *r_E, *r_H, *g_E, *g_H,integer_bit;
  18.     double length, x, y, lamda_design, lamda, tmp;
  19.     struct FCOMPLEX temp_E, temp_H, e_gain, h_gain;
  20.     /* need to allocate space for FCOMPLEX types. Since there is no Numerical 
  21.     Recipes routine, I'll use the standard method. Since I've always used
  22.     elements positions 1 to N, I'll just wast the extra location */
  23.     r_E=dvector(1L,(long) elements); 
  24.     r_H=dvector(1L,(long) elements);
  25.     g_E=dvector(1L,(long) elements);
  26.     g_H=dvector(1L,(long) elements);
  27.     e_gain.r=0;   /* set to zero. Necessary as we sum into these */
  28.     e_gain.i=0;
  29.     h_gain.r=0;
  30.     h_gain.i=0;
  31.     /* convert theta and thi to radians from degrees */
  32.     theta=theta*M_PI/180;
  33.     phi=phi*M_PI/180;
  34.     lamda_design=3e8/design_frequency;
  35.     lamda=3e8/actual_frequency;
  36.     for(i=1;i<=elements;++i) 
  37.     {
  38.         length=coordinates[i].length/lamda;
  39.         x=coordinates[i].x/lamda_design;
  40.         y=coordinates[i].y/lamda_design;
  41.  
  42.         /* for E -plane */
  43.         r_E[i]=sin(theta)*x;
  44.         if(fabs(theta) < TINY)    /* avoid division by zero if theta=0 */
  45.             g_E[i]=0;
  46.         else
  47.             g_E[i]=(cos(M_PI*length*cos(theta))-cos(M_PI*length))/sin(theta);
  48.  
  49.         /* for H -plane */
  50.         r_H[i]=x*cos(phi)+y*sin(phi);
  51.         g_H[i]=1-cos(M_PI*length);
  52.         /* printf("g_H[%d]=%.16lf \n", i, g_H[i]);         */
  53.         temp_E.r=0;
  54.         temp_E.i=2*M_PI*r_E[i]*F;
  55.         temp_E=E_to_complex_power(temp_E); /* exp(j 2 pi r[i] F) */
  56.  
  57.         temp_H.r=0;
  58.         temp_H.i=2*M_PI*r_H[i]*F;
  59.         temp_H=E_to_complex_power(temp_H);
  60.       /* printf("element %d temp_H.r=%lf temp_H.i=%lf\n",i, temp_H.r, temp_H.i); */
  61.         /* get element currents */
  62.         temp_E=Cmul(temp_E,current[i]);
  63.         temp_H=Cmul(temp_H,current[i]);
  64.         /* printf("element %d: current = %.10lf i%.10lf = %.10lf (mag) at %lf degrees\n", i,current[i].r, current[i].i, sqrt(current[i].r*current[i].r+current[i].i*current[i].i), 180*atan2(current[i].i,current[i].r)/M_PI);         */
  65.       /* printf("element %d I temp_H.r=%.10lf temp_H.i=%.10lf\n",i, temp_H.r, temp_H.i); */
  66.         e_gain=Cadd(e_gain,RCmul(g_E[i],temp_E));
  67.         h_gain=Cadd(h_gain,RCmul(g_H[i],temp_H));
  68.         /* printf("h_gain=%.16lf + i %.16lf \n", h_gain.r, h_gain.i); */
  69.     }
  70.     /* there will be a divide  by zero in calculating
  71.     equ 1.28, of Lawsons book, (see page 1-12)
  72.     at theta=0,180,360 degree, infact anywhere where 
  73.     sin(theta) is zero */
  74.  
  75.     if( modf(theta/M_PI,&integer_bit) < TINY || fabs(theta-M_PI)<TINY)
  76.     {
  77.         *gain_E_plane=-999;   /* a very very small number */
  78.     }
  79.     else
  80.     {
  81.         tmp=(Cabs(e_gain)*Cabs(e_gain)*60/pin);  
  82.         *gain_E_plane=10*log10(tmp);  
  83. #ifdef DEBUG
  84.     if(errno)
  85.     {
  86.         fprintf(stderr,"Errno =%d in gain() of gain.c after E\n", errno);
  87.         printf("tmp=%lf pin=%lf theta=%lf gainE=%lf\n", tmp, pin,theta,*gain_E_plane);
  88.         exit(1);
  89.     }
  90. #endif
  91.     }
  92.     tmp=(Cabs(h_gain)*Cabs(h_gain)*60/pin);  
  93.     *gain_H_plane=10*log10(tmp);  
  94.  
  95. #ifdef DEBUG
  96.     if(errno)
  97.     {
  98.         fprintf(stderr,"Errno =%d in gain() of gain.c after H\n", errno);
  99.         printf("tmp=%lf pin=%lf \n", tmp, pin);
  100.         exit(1);
  101.     }
  102. #endif
  103.     free_dvector(r_E,1L,(long) elements); 
  104.     free_dvector(r_H,1L,(long) elements); 
  105.     free_dvector(g_E,1L,(long) elements); 
  106.     free_dvector(g_H,1L,(long) elements); 
  107.  
  108. #ifdef DEBUG
  109.     if(errno)
  110.     {
  111.         fprintf(stderr,"Errno =%d in gain() of gain.c\n", errno);
  112.         exit(1);
  113.     }
  114. #endif
  115. }
  116.  
  117. struct FCOMPLEX E_to_complex_power(struct FCOMPLEX x)
  118. {
  119.     struct FCOMPLEX a, answer;
  120.     double real;
  121.  
  122.     /* First real part */
  123.     real=exp(x.r);
  124.     /* now the imaginary part */
  125.     /* Exp(i theta) = cos(theta) + i sin(theta) */
  126.     a.r=cos(x.i);
  127.     a.i=sin(x.i);
  128.     answer=RCmul(real,a);
  129.     return(answer);
  130. }
  131. #undef TINY
  132.